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We develop a 3 dimensional computer code to study a coalescing neutron star binary. 
The code can currently follow the evolution up to two stars begin to merge from two spherical 
O ■ stars of mass IMq and radius 8.9km with separation 35.4km. As for coordinate conditions, 

we use conformal slicing and pseudo-minimal distortion conditions. The evolution equations 
' for the metric is integrated using the CIP method while the van Leer's scheme is used to 

integrate the equations for the matter. We present a few results of our simulations including 
gravitational radiation. 



^ I §1. Introduction 

00 

O ■ We report on the present status of our computer code development to study the 

fully general relativistic evolution of spacetimes and matter. We solve the complete 
set of the Einstein equations in the absence of any symmetries. The main motivation 
Q\ ■ of this project is studying promising sources of strong gravitational waves, such as 

coalescing neutron star or black hole binaries. The first detection by an interfero- 
metric detector of gravitational waves such as LIGO, VIRGO, GEO600 and TAMA 
q may be the waves from a black hole binary (see the contributions of K.S. Thorne, 

S. Kawamura and B. Schutz in this volume), but the astrophysical importance of a 
coalescing neutron star binary is not smaller than a black hole binary. A detailed 
analysis of the detected waves near the merger of two stars, for example, will give 
information on the size of a neutron star and then on the equation of state of high 
density matter. On problems other than the gravitational wave physics, study on 
process of coalescence of a neutron star binary is helpful for models of a source of a 
gamma ray burst (see the contribution of P.I. Meszaros in this volume). 

A compact star binary in the earliest stage is almost stationary but gravitational 
radiation reaction makes it evolve in the adiabatic manner. When the separation 
between stars is large enough, the stars can be regarded as point particles. As the 
separation is getting small, however, the size effect of the stars and the relativistic 
effects become important. Thus the complete set of the Einstein equations and 
the general relativistic hydrodynamics equations for a neutron star binary must be 
solved. In addition, the process of coalescence of a binary is not axisymmetric. Thus 
a three dimensional (in space) code is necessary and it requires a great power of 
computers. In 1994, we became able to use a supercomputer with memories of tens of 
GBytes and a peak speed of nearly a hundred GFLOPS. With this computer, we can 
manage a 200 3 grid. Then we started to make a 3D, fully general relativistic code for 
a coalescing neutron star binary. In this paper, we show the present status of our code 



typeset using ■pTpTfX.sty <ver.0.8> 



2 



development. The present code does not include a lot of realistic physics processes 
such as neutrino emission and a realistic equation of state near and beyond the 
nuclear density. So it is a milestone of 3D numerical relativity. In §2, we describe the 
mathematical framework of the calculation, including constraint equations, evolution 
equations of the metric, relativistic hydrodynamics equations, coordinate conditions 
and evaluation of the gravitational waves. In §3, numerical methods we use for 
elliptic equations and evolution equations are shown. In §4, we present numerical 
results of simulations of coalescing neutron star binaries obtained using the present 
code. In §5, we summarize and outline future work planned for this code. 

§2. The Mathematical Framework 

2.1. Constraint equations 

In the (3+1) (or ADM) formalism of general relativity, the spacetime metric is 
written as 

ds 1 = -a 2 dt 2 + 7i J -(dx i + p i dt)(dx j + ftdt), (2-1) 

where a, (3 l and 7^ are the lapse function, the shift vector and the intrinsic metric 
of three-metric. The Einstein equation reduces to the four constraint equations 

R + K 2 - KijK ij = 16tt Ph , (2-2) 

Dj (K j i - 5 j iK^j = 8vr J t (2-3) 

and 12 evolution equations 

dan = -2aK i:j + DiPj + Djfii, (2-4) 



dtKij = 



DiDja 



Rij - 8vr {Stf + i 7ij (j) H - S e e ) } 
+q (KKij - 2K a K t j ) + K et Djf3 e + K ti D^ + ^D t K^ (2-5) 



where Rij is the Ricci tensor, Kij the extrinsic curvature tensor, K = K*t, and Di the 
covariant derivative associated with 7^ . The quantities p H , Ji and Sij are the energy 
density, the momentum density and the stress tensor, respectively, measured by the 
observer moving along the line normal to the spacelike hypersurface of t =constant. 

In order to give an initial data, we should find a three-metric and extrinsic 
curvature which satisfy the constraint equations (2-2) and (2-3). It means that the 
constraint equations are solved giving p H and Jj. Here we assume that K = and 
jij is conformally flat at t = 

lij = <t> 4 7ij, (2-6) 



where 7^ is the flat space metric. Defining the conformal transformation as 

(2-7) 



p B = 4> % p H and Ji = (p 6 Ji, 
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Eq.(2-3) reduces to 

DjK'i = SirJi, (2-8) 

where D{ is the covariant derivative associated with jij. 1 ^ The traceless extrinsic 

curvature can be decomposed with the transverse traceless part K^ T and the longi- 
tudinal traceless part {LW)^ 2 ^; 

K ij = K? T + (LW) ij , (2-9) 

where 

(LW)ij = D t Wj + DjWi - | 7ij&W t . (2-10) 
Assuming K^ T = 0, Eq.(2-8) becomes 



AWj + ^DiD j Wj = SirJi, (2-11) 



where A = D % D; L . Equation (2-11) is the coupled elliptic equation but can be reduced 
to four decoupled Poisson equations: 

A X = 6ttcW, (2-12) 

AWi = SnJi - ^d iX , (2-13) 

where x = diW{. 

The Hamiltonian constraint equation (2-2) reduces to the nonlinear Poisson 
equation 

Lcj) = - l<f>~ 7 KijK ij . (2-14) 

2.2. Evolution equations of the metric 
Defining the following variables 

4> = (det( 7ij ))^ , lij = <t>~Sij , l ik lkj = (2-15) 

Kij = STF^Kij), K*j = rf k K kj , K l i = ^ k K\, (2-16) 

where 

STF(Ai) = \ (A l3 + Aji - lw H A k ^j , (2-17) 

the evolution equation (2-5) reduces to 

dtKij - pdeKij = <p- 4 [STF {a (R tJ - 8irS tJ ) - D t D ja }} 

+ a(KK ij -2K ii K e j ') (2-18) 



d t K - p e d e K = a 



-DWia, (2-19) 
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and Eq.(2-4) to 



dtjij = -2 



aK ij -STF{D i (<!>'%)} 



— Aij 



or 



dtlij - P e d e jij = -2 
The conformal factor <\> obeys 



aKij - STF (judjP*) 



or 



d t <t>-p e d t <t> = -£(aK-dtP e ) 

/ 5 



(2-20) 
(2-21) 

(2-22) 
(2-23) 



Equation (2-22) follows from the trace of Eq.(2-4), while Eq.(2-23) from the Hamil- 
tonian constraint (2-2). 



2.3. Relativistic hydrodynamics equations 

We assume the perfect fluid stress-energy tensor, which is given by 

= (p + pe + p)u fM u v + pg^ u , 



(2-24) 



where p, e and p are the proper mass density, the specific internal energy and the 
pressure, respectively, and is the four-velocity of the fluid. The energy density 
p H , the momentum density Jj and the stress tensor Sij of the matter measured by 
the normal line observer are, respectively, given by 



p H = nWT^, Ji = -hfrfT^ and 5^- = hfhfT^, 



(2-25) 



where is the unit timelike four-vector normal to the spacelike hypersurface and 
h^ u is the projection tensor into the hypersurface defined by 



(2-26) 



The relativistic hydrodynamics equations are obtained from the conservation of 
baryon number, V jU (pu /1 ), and the energy-momentum conservation law, V ' V T^ ' . In 
order to obtain equations similar to the Newtonian hydrodynamics equations, we 
define p N and uf as 



p N = ^fyau°p and u 



Ji 



n ' 
au u p 



(2-27) 



respectively, where 7 = det^). Then the equation for the conservation of baryon 
number takes the form 

d t p N + d e (p N V e ) =0, (2-28) 

where 



u 



V t = ^ = 



P + Pl 



-a 1 . 



(2-29) 
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The equation for momentum conservation is 

dt{p N uf) + d e (p N ufV^ = -yfjadip - ^{p + p H )dia 

The equation for internal energy conservation becomes 

d t(p N e) + d e {p N eV^j = -pd u (^au v ) . 
To complete hydrodynamics equations, we need an equation of state, 

P = p{e,p)- 



(2-30) 



(2-31) 



(2-32) 



The right-hand side of Eq.(2-31) includes the time derivative. For a polytropic 
equation of state, p = (T — l)pe, however, the equation reduces to 



dt(p N e N ) + d e (p N e N V £ ) = -p N d e V e , 



where 



r-i 



p N = (r- i)p n £ n = (Viau ) P- 



(2-33) 

(2-34) 
(2-35) 



2.4. Coordinate conditions and gravitational waves 

The choice of a shift vector (3 l and a lapse function a is important because the 
stability of the code largely depends on it and because it is intimately related to 
the extraction of physically relevant information, including gravitational radiation, 
in numerical relativity. 

Since the right-hand side of Eq.(2-20), Ay is trace-free, the determinant of 7y 
is preserved in time and the condition 



D J A i:j = 



(2-36) 



produces the minimal distortion shift vector. 3 ) It must be a good choice of the spatial 
coordinate, but it is too complicated to be solved. Instead we replace the covariant 
derivative by the partial derivative 



d{Aij — 0, 

which yields an elliptic equation for the shift vector (3 l : 



V 2 /f + x -did^ 



2d, 



aK l3 - STFihudjF) - -tfdihij 



where 



h 



(2-37) 

(2-38) 
(2-39) 
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We call this condition as the pseudo-minimal distortion condition. 

As for slicing condition, we choose the conformal slicing, 4 ) where the lapse func- 
tion a is given by 



with 



a = exp 
1. Since in this slicing, 



1 



a 



1 + . 



i _ M 

2r_ 

1 + M 

1 ^ 2r 



(2-40) 



(2-41) 



for large r, the space outside the central matter quickly approaches the Schwarzschild 
metric in the isotropic coordinates. Thus hij defined by Eq.(2-39) can be considered 
as gravitational wave parts for large r and the total energy of the gravitational waves 
is given by 



=I dt L 



dQ r 2 Pi 



r— >oo 
1 



GW 



^gw ~~ 2,2TT AijAij ' 



(2-42) 
(2-43) 



since Eq.(2-37) means that Aij is the transverse-traceless part of the time derivative 
of 7jj for large r. 



§3. Numerical Methods 



3.1. Coordinate system 

As for the spatial coordinates, we use a Cartesian coordinate system (x,y, z). 
We solve partial differential equations (PDEs) using a finite-difference scheme, where 
a finite grid is introduced and the PDEs are replaced by finite-difference approxi- 
mations. It may be convenient to use a spherical coordinate system (r, 9, (f) with 
regard to saving of computer memory and an outer boundary condition at large dis- 
tance from a compact source. But the use of spherical coordinates leads to a vexing 
problem of coordinate singularities at the origin and along the polar axis. 

3.2. Elliptic equations 

The constraint equations reduce four Poisson equations (2-12) and (2-13), and 
a nonlinear Poisson equation (2-14). These equations are solved once for initial 
data. The shift vector is given by the coupled elliptic equation (2-38), which is 
solved at every time step. The conformal factor <f> is, in turn, determined either by 
the hyperbolic equation (2-22) or by the elliptic equation (Hamiltonian constraint) 
(2-23). However, since it is known that the conformal slicing yields instability in <f> 
and a, 5 ) we solve the Hamiltonian constraint (2-23) at every time step to prevent 
the onset of the instability caused by numerical errors from occurring quickly. Thus 
the solution of the elliptic equations consumes the greatest part of the CPU time. 

A non-linear equation such as 



A0(t) = -47rS(«Kt)), 



(3T) 
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where S(cp) is a non-linear function of (f>, must be solved by an iterative method: 

= A -i s ^(i)) for 7 = 0, 1,2,- ••. (3-2) 

To solve Eq.(2T4) for initial data, we should made the iteration. For Eq.(2-23), 
however, we don't made the iteration to reduce the required CPU time. Instead, the 
right-hand side S((f>(t)) is calculated using an approximate value of </>(i) evaluated by 
extrapolation from the values of <fi at the previous two time steps t — At and t — 2 At. 
We examined the difference of the solution using the extrapolated source from the 
solution obtained by the iterative method, and found that it is small enough if At 
is as small as we use for our actual calculations. 

We use the extrapolated value of [3 l in calculating the right-hand side of Eq.(2-38), 
too. Note that terms including (3 l vanish and we don't need the extrapolated value 
at t = 0, since hij = 0. Thus Eq.(2-38) reduces four Poisson equations (2-12) and 
(2-13), replacing W, to (3\ 

Now all the elliptic equations reduce Poisson equations, which are solved using 
a pre-conditioned conjugate gradient method. 6 ) The boundary condition of <f> giving 
by Eq.(3-1) is that 

, M d x x k „ ( 1 \ . . 

for large r, where 

M = Jsd 3 x, d k = Jsx k d 3 x. (3-4) 
The boundary condition of \ an d Wi (or f3 l ) giving by Eqs.(2-12) and (2-13) is that 6 ) 
PiX* 3M U VM^xi n fl\ , , 

w = Pkxix ' _ m 

(7Mjj - Mjj - M kk 5jj) xi 3%^iV |Q /1\ (3g) 



for large r, where 

Pi = j Jid 3 x, Mij = J Jix j d 3 x. (3-7) 

3.3. Evolution of matter 

Hydrodynamics equations Eqs.(2-28), (2-30) and (2-33) are written as 

d t Q + dz(QV e ) = F, (3-8) 

which is split into two phases: 

dtQ + d t (QV e ) = 0, (3-9) 



d t Q = F. 



(3-10) 
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Here we call Eq.(3-9) the advection phase and Eq.(3-10) the non-advection phase. 
The advection phase is solved using van Leer's scheme 7 ) to obtain Q^j^* from Q^ bc 
and (V e )™ bc , where Q™ bc is the value of Q at n-th time step t = t n and a spatial grid 
point (x a ,y b , z c ). If the non-advection phase is solved as 

QZi 1 = Qijk 1] * + At • F ijk (3-n) 

with Ff- k calculated using quantities at t = t n , it is of first-order accuracy in time. 
To achieve second-order accuracy, we use (1) a two-step algorithm 

Qtk h = Q { S i> + \ At - F ^ ( 3 - 12 ) 

= Q { St X> + (3-13) 

n+— u-\- — 

where Q^ k 2 is used only in calculating F^ k 2 , or (2) a extrapolated source algorithm 

yn+l _ n (n+l)* , A , 



where F™^ 2 given by 



F^=FZ k + \[F% k -F^) (3-15) 

is the value of F at t = t n + \At evaluated by extrapolation from the values at t = t n 
and t = t n — At. Since some test calculations indicated that these algorithms give 
almost the same results if At is small, we will use the extrapolated source algorithm. 

Since the outer boundary is sufficiently far from the region of the matter distri- 
bution, we can safely set Q = at the outer boundary. 

3.4. Evolution of metric 

The evolution equations for Kij (2-18), K (2-19) and jij (2-21) can be written 

as 

d t Q + v e d £ Q = F, (3-16) 

with v e = —j3 e . It is solved using the CIP (Cubic-Interpolated Pseudoparticle/Prop- 
agation) method. 8 ) In the CIP method, both Eq.(3T6) and its spatial derivatives 

dt(d a Q) + v%(d a Q) = -{d iV e )(d e Q) + 8 a F (3-17) 

are solved. Numerical diffusion during propagation of Q can be reduced, since the 
time evolution of Q and its derivatives are traced. Equations (3-16) and (3-17) are 
split into the advection and the non-advection phases similarly to Eqs.(3-9) and 
(3-10). In the non-advection phase, F and d a F at t = t n + \At are evaluated by 
extrapolation from the values at previous two time steps. 

By the way, special care must be taken to calculate the Ricci tensor appearing 
in the right-hand side of Eq.(2T8). With conformal transformation Eq.(2T5), the 
Ricci tensor Rij associated with 7^ is given by 

Rij = + Rif\ (3- 18) 
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where 



+24>- 2 (siPifiiDjt) - %(D k <l>)(D k <i>)) 



and Rij is the Ricci tensor associated with 7^, which is given by 
1 



Rij 



1 M {hi,jk + hej,ik - hij,kt) + V^.fc (%J + hej,i - hij 



ij,e) 



k 



J-\K 7-1. 

i til , 



jk- 



(3-19) 



(3-20) 



The second derivatives of are replaced by finite differences in numerical calcula- 
tion but the numerical precision of terms such as hij^e with k / £ is not so good, 
while the degree of precision in calculation of hij^e with k = I is the same as that 
of the first derivatives. Inaccuracies in hij ki will cause a numerical instability. The 
'pseudo-minimal distortion condition', Eq.(2-37), however, leads S^h^^ = and 
therefore 

J k (hgijk + h£j t ik - hij t ki) = -hij^k + tp he {hujk + Wj,ik - hij t ki) , (3-21) 

where tp^ = 7 iJI — 8 % K Inaccuracies in numerical values of ip k ^h-£i,jk wu l n °t affect 
integration of Eq.(2-18) so seriously since both tp and hgi jk are smaller than unity. 

Boundary conditions for Kij and 7^ are important, because the boundary is 
not so far from the stars and therefore inappropriate boundary conditions cause 
reflections of the outgoing waves. We apply outgoing conditions for Ky and 7^, 
since non-wave parts decrease rapidly in our coordinate conditions. The condition 
can be written as 

H(at - <j) 2 r) 



Q{t,x l 



where r = \Jx 2 + y 2 + z 2 . Here H satisfies an advection equation 

8 t H + c e d e H = 0, 

where 



ax 



c = 



h2 • 



It implies 



d t Q + c e d £ Q 



a 

rcf> 2 



(3-22) 

(3-23) 
(3-24) 
(3-25) 



Equation (3-25) is solved along with Eq.(3-16) using the CIP method. 



§4. Numerical Simulations of Coalescing Binary Neutron Stars 

Now we show numerical results on simulations of coalescing binary neutron stars. 
As for a equation of state, we use the 7 = 2 polytropic equation of state. To 
reduce the required memory, we assume the symmetry with respect to the equatorial 
(z = 0) plane. A uniform grid of size 191 x 191 x 96 with Ax = Ay = Az = 
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Fig. 1. Initial configuration of a coalescing neutron star binary. 

lGM Q /c 2 = 1.5km is used. The numerical boundary is located at 95GM & /c 2 from 
the origin. In order to keep precision, we set the time step At as 0.01GM Q /c 3 = 
5 x 10~ 8 seconds. Calculations are performed on Fujitsu VPP500/80 at High Energy 
Accelerator Research Organization (KEK) and Fujitsu VPP300/16R at the National 
Astronomical Observatory (NAO), Japan. The required memory is approximately 
8GBytes. 

4.1. Initial data 

We performed calculations with two kinds of initial data. One is the corotation 
(synchronized rotation) case as for ordinary stellar binaries, in which the momentum 
density Jj appearing in Eq.(2-8) of each star is given by 

Ji = P B e ijk f2 j x k , (4-1) 

where is the complete antisymmetric symbol and the angular velocity W is 
constant. The other is the irrotational (zero-circulation) for gravitational 

radiation driven compact stellar binaries, in which Ji is given by 

Ji = P B e ijk^ j x k c , (4-2) 

where x\ is the position of each star and is constant for an individual star. In both 
cases, we assume that the binary consists of identical spherical stars of rest mass 
1.0M Q , radius 8.9km located at (x c ,y c ,z c ) = (±17. 7km, 0, 0). The orbital angular 
momentum Q points to the z-direction. 

4.2. Corotation case 

The parameters of initial data for the corotation case are shown in the second 
column of Table I. It take about 10 hours on VPP500 using 64 processor units up 
to 26000 time steps. Figure 2 shows the evolution of the density and the velocity on 
the x-y plane. Movies for it including the evolution on x-z and y-z planes can be 
obtained from our Web site. 9 ) At the ~26000th time step, the onset of instability 
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Table I. Parameters of initial data. The angular momentum parameter q is Jt/(GM^ oi /c). 





Corotation 


Irrotation 


Orbital angular velocity O [sec -1 ] 


2.4 x 10 :l 


2.4 x 10 ;i 


Angular momentum J t [GMq /c] 


3.9 


3.5 


Collision velocity v x [c] 


0.074 


0.071 


Total gravitational mass M^ i [Mq] 


1.96 


1.95 


Angular momentum parameter q 


1.03 


0.93 



on <fi in the central region occurs and it grows exponentially in a few time steps. At 
the time, only 3 grid points or so are included from the center to the surface of the 
star. So that numerical errors grow on account of unstable modes resident in the 
conformal slicing. 

A thick disk around the merged star is likely to be formed as shown in Fig. 3. 
The mass of the disk is about a few percent of the total mass. In addition, there 
are indications of outflow along the rotation axis. However, since the precision of 
the numerical calculation at this time is not so good, these facts must be confirmed 
with more precise calculation. 

The propagation of the "energy density of gravitational waves" r 2 p GW is shown 
in Fig. 4, movies of which are also provided on our Web site. 9 -* A spiral pattern 
appears on the x-y plane, while different patterns with peaks around z-axis appear 
on the x-z plane. This can be explained naively by the quadrupole wave pattern 
given by 

r 2 

r 2 p = (Aijf oc cos 2 + sin 2 9 sin 2 (2f2 (t - r) - 2<p)/4. (4-3) 

327T 

On the x-y plane, where 9 = tt/2, p GW is constant along the spiral of r + <p/£2 = 
constant, while near z-axis, where 9 « 0, P GW oc cos 2 9. 

The waves pass through the numerical boundary with little reflection. It means 
that the boundary condition Eq.(3-25) behaves well. However, a small reflection on 
the x-, y- and z-axes interfere with the outgoing waves and interference patterns 
appears after t ~ 0.8msec. 

Figures 5 shows the waveforms on the z-axis monitored between 75M and 
85Mq from the origin. The emitted waveform falls off with r _1 as expected. The 
luminosity evaluated between 75M & and 85M as functions of t and the total energy 
of the gravitational waves up to t are shown in Fig. 6. Comparing Figs. 4, 5 and 
6 with Fig. 2, the main part of the waves begins to be emitted when the merger 
develops to a certain extend, but only the waves up to the time when the merger 
begins have arrived at the numerical boundary. The minor peak around t — r = 
corresponds the waves included unexpectedly in the initial data. 

The amplitude of gravitational waves is 

h = 5x 1(T 22 ( ^ ) (4-4) 

VIOMpcy v ' 

and the total energy emitted is 

E GW (t < 1msec) = l(T 3 M c 2 = 2 x 10 51 erg ~ 0.05% of M tot c 2 . (4-5) 






. 2. Density p N and velocity V 1 on the x-y plane for the corotation case. Time in units of 
milliseconds is shown. Solid lines are drawn from p m j n through pmax in steps of p m ; n , where 
Pmin = 8 x 10 14 g-cm~ 3 and pmax = 4 x 10 16 g-cm -3 . Dashed lines are at 1/2, 1/4, 1/8, 1/16, 
1/32 of p min . 
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. 3. Density p N and velocity V % in the central region on the x-y, x-z and z-y planes for the coro- 
tation case. Parameters for contour lines are the same as Fig. 2 but pmax 
and p min = O.lpmax- 



3.41 x 10 16 g-cm~ 3 
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Fig. 4. The propagation of the gravitational waves for the corotation case. The "energy densities 
of gravitational waves" r 2 pQ-yy on the x-y and x-z planes are shown as gray scale figures up to 
the numerical boundary. Time in units of milliseconds is shown. The upper figure of each time 
is on x-z plane and the lower one is on x-y plane. 



along the z-axis 




54 

X 10 erg/sec 



0.5 
t - r (msec) 

Fig. 5. Waveforms h+ and h x observed 
on the z-axis for the corotation case. 
They have been monitored at r = 
75M ,8OM and 85M Q . 




0.0 0.5 

t - r (msec) 



Fig. 6. Luminosity (in units of erg/sec) and to- 
tal energy (in units of erg) of the emitted 
gravitational waves. They are evaluated at 
r = 75M , 8OM and 85M . 
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Fig. 7. Density p N and velocity I/ 1 on the x-y plane for the irrotation case. Parameters are the 
same as Fig 2. 



The angular momentum of the emitted gravitational waves is given by 

dJ 



AJ = 



dt 



dt, 



to be 



Z\J(t < 1msec) = 1.8 x W' z GM^/c ~ 0.5% of J t . 
4.3. Irrotation case 



(4-6) 
(4-7) 
(4-8) 



The parameters of initial data for the irrotation case are shown in the third 
column of Table I. The evolution of the density and the velocity on the x-y plane 
are shown in Fig. 7. 9 ) 

The propagation of the "energy density of gravitational waves" r 2 /? GW are shown 
in Fig. 9, 9 ) the waveforms on the z-axis in Fig. 10 and the energy of the waves in 
Fig. 11. 

The evolution is in principle the same as the corotation case, but no outflow 
along the rotation axis appears in Fig. 8. It is because the merger is progressing 
slower. It is the case in post-Newtonian calculations also that the progress of the 
merger for the irrotation case is slower than for the corotation case. 6 ) It causes the 
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Fig. 8. Density p N and velocity V" in the central region on the x-y, x-z and z-y planes for the 
irrotation case. Parameters are the same as Fig. 3. 
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Fig. 9. The propagation of the gravitational waves for the irrotation case. The "energy densities 
of gravitational waves" r 2 on the x-y and x-z planes are shown similarly to Fig. 4 
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Fig. 11. Luminosity (in units of erg/sec) and to- 
tal energy (in units of erg) of the emitted 
gravitational waves. They are evaluated at 
r = 75M Q , 8OM and 85M Q . 



emission of gravitational waves to decrease somewhat. The amplitude, the energy 
and the angular momentum of the gravitational waves emitted up to 1msec are, 
however, almost the same as in the corotation case: 

h = 5 x 1(T 22 ( -]—) , 

VioMpcy 

E GW (t < 1msec) ~ 0.05% of M tot c 2 , 
AJ{t < 1msec) ~ 0.5% of J t . 

§5. Summary 

We have described a numerical code for simulations of coalescing binary neutron 
stars and presented some results with it. This code was developed after extensive 
study of a post-Newtonian code for coalescing binary neutron stars. 6 ) To solve hy- 
drodynamics equations and elliptic equations, the general relativistic code uses the 
same methods as the post-Newtonian code does. For evolution equations for the met- 
ric, however, we adopted the CIP method, which made the evolution of the metric 
stable even with a coarse grid. 

We faced instability when coalescence of two stars advanced to be probably a 
black hole. It is caused by a characteristics of the conformal slicing as well as by the 
coarseness of the grid. This fact may suggest that other slicing such as the maximal 
slicing must be adopted. Thus we are now making another code with the maximal 
slicing. With the maximal slicing, we must solve a elliptic equation for a, which 
is obtained by setteing the right-hand side of Eq.(2T9) equal to zero. The cost of 
solving this equation is great. In addition, to catch the wave form of gravitational 
waves, we have to apply a gauge-invariant wave extraction technique. 10 ) Since new 
super computer, 7 times faster with 20 times greater memory, will be available 



(4-9) 
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at KEK from the beginning of 2000, more precise calculation with the maximal 
slicing will be carried out. To calculate the right-hand side of Eq.(2-38), we used the 
extrapolated values for (3 l as describe in §3.2, but the precision of the extrapolation 
will reduce in the highly dynamical phase of coalescence. We will be able to solve 
Eq.(2-38) more precisely by an iterative method on the new machine. 
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